Stability of dark solitons in a Bose-Einstein condensate trapped in an optical lattice 
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We investigate the stability of dark solitons (DSs) in an effectively one-dimensional Bose-Einstein 
condensate in the presence of the magnetic parabolic trap and an optical lattice (OL). The analysis 
is based on both the full Gross-Pitaevskii equation and its tight-binding approximation counterpart 
(discrete nonlinear Schrodinger equation). We find that DSs are subject to weak instabilities with 
an onset of instability mainly governed by the period and amplitude of the OL. The instability, if 
present, sets in at large times and it is characterized by quasi-periodic oscillations of the DS about 
the minimum of the parabolic trap. 



The experimental creation and great advancement in 
the theoretical understanding of Bose-Einstein conden- 
sates (BECs) [1] have stimulated a lot of interest in non- 
linear matter waves, including dark [2] and bright [3] soli- 
tons. The dynamics of dark solitons (DSs) in the pres- 
ence of the external magnetic trap has been extensively 
studied [4], including thermal [5] and dynamical [6] in- 
stabilities. More recently, apart from the rectilinear DSs, 
ring-shaped counterparts were predicted in BECs [7] . 

The study of nonlinear excitations is particularly rel- 
evant for BECs trapped in optical lattices (OLs) gener- 
ated by interference patterns from laser beams illuminat- 
ing the condensate [8-13]. The controllable character of 
the OL allows for the observation of numerous phenom- 
ena, such as Bloch oscillations [10,14] and Zener tunnel- 
ing [8] (in the presence of an additional linear external 
potential), or classical [15] and quantum [13] superfluid- 
insulator transitions. 

Apart from the mean-field description via the Gross- 
Pitaevskii (GP) equation [8-13], a BEC trapped in a 
strong OL may be described, in the tight-binding limit, 
by the discrete nonlinear Schrodinger (DNLS) equation 
[16]. This approximation is not always accurate, but 
its applicability can be systematically examined [17]. In 
cases where such a reduction is possible (e.g., when the 
chemical potential is much lower than the height of the 
potential barriers induced by the OL), the DNLS model 
is particularly relevant and has been successfully applied 
in many instances (for a recent review on DNLS see, e.g., 
[18] and references therein). 

In this paper, we study DSs in repulsive BECs (i.e., 
positive-scattering-length collisions) in the presence of 
OLs. We use both the continuous-GP and DNLS equa- 
tions. In particular, we assume a cigar-shaped BEC, 
which can be described by the following normalized 
quasi-one-dimensional GP equation [1,19,20], 

iut = —Uxx + |u|' u -|- \kx + Vo sin (27ra;/A)] u . (1) 

Here, u{x,t) is the mean-field wave function, while the 
terms in the square brackets represent the external mag- 



netic trap and the OL potential, respectively, with the 
strengths k and Voj while A is the wavelength of the in- 
terference pattern created by the laser beams. 

To study the dynamics of a DS in the framework of Eq. 
(1), we consider an initial condition similar to an ansatz 
proposed for the description of DSs in BECs in Ref. [21] 



uq{x) ~ u^y{x) tanh(a; — xq), 



(2) 



where utf ~ ■\/max(0, /i — kx'^) is the Thomas-Fermi 
(TF) expression for the background wave-function dis- 
tribution [1] (/x is the chemical potential) and xq is the 
initial location of the DS's center. In most cases, we set 
xo = 0, i.e., the dark soliton is placed at the bottom of 
the magnetic trap. 

In the tight-binding limit, Eq. (1) reduces to the fol- 
lowing DNLS equation [16], 



-C(u„+i +' 



- 2u„) + |u„|^u„ -I- kn? 



(3) 



where the dot denotes time derivative, n is the lattice- 
site index, and C is the so-called coupling constant (see 
e.g., Refs. [16,17] for exact expressions and relevant esti- 
mates). An initial condition for a DS in the case of Eq. 
(3) can be given by a straightforward discretization of the 
continuum ansatz (2). Typically, simulations were run 
for a lattice with 200 sites, and free boundary conditions 
were used for both the continuum and discrete models. 
In fact, it has been verified that the results are insensitive 
to the choice of boundary conditions. Note that DSs in 
discrete lattices (in the absence of a parabolic trap) were 
already studied in the framework of the DNLS equation 
[22], revealing that they are subject to oscillatory insta- 
biUties [23]. 

Stationary solutions to Eq. (3) are sought for in the 
form M„ = exp(— i^i)w„, where fx is the chemical poten- 
tial [1], which leads to the steady-state equation, 

-C{Vn+l + Vn-1 - 2Vn) + (|u„|^ -l- kr? - ^)u„ = 0, (4) 

for a discrete function w„ . Once a solution of Eq. (4) is 
found, linear stability analysis is performed by looking 
for perturbed solutions of the form 



-ifj,t 



■eanC 



iut 



ebr, 



-iLj*t'\ 



(5) 



where the asterisk denotes complex conjugation. The 
ensuing eigenvalue problem, for the eigenfrequency- 
eigenfunction pair {a;, {a„, 5*}}), is then solved numeri- 
cally with uj ~ ujr + i^i (where the subscripts denote the 
real and imaginary parts). In what follows, since Eqs. 
(1) and (3) admit an additional rescaling, the chemical 
potential has been fixed at /i = 1. 




FIG. 1. The left and right top panels show, respectively, 
the (square root of the) number of atoms in the condensate 
N and the largest unstable eigenfrequency uii for the discrete 
dark soliton as a function of k for C = 1. The profile of 
the (exact) solution is shown, for k — 0.001 and k — 0.119 
(just prior to the termination of the branch), in the left mid- 
dle and bottom panels, respectively. The right middle and 
bottom panels show the spectral plane {u>r,i^i) of the corre- 
sponding eigenfrequencies, the subscripts standing for the real 
and imaginary parts. 

We first consider the DNLS model with a fixed cou- 
pling constant C — 1 while the strength k of the parabolic 
trap is varied. In this case, it is found that the assumed 
configuration, in the form of a DS on top of the TF back- 
ground, exists only up to a critical value of fccr ~ 0.12, 
as is seen from the top left panel of Fig. 1, which com- 
plies with the well-known fact that the DS cannot have a 
width essentially smaller than the healing length [1], and 
thus cannot exist in a very narrow trap. 

This solution family is quantified by the top panels of 



Fig. 1. The norm N — \\v„ 



Sn ^n (number of atoms 



N) as a function of k is depicted in the top left panel. 
The dependence of the largest instability growth rate uji 
on fc, which is shown, for the same discrete-DS solution, 
in the top right panel of Fig. 1, reveals a narrow stability 
window close to fe = (for 0.002 < A; < 0.004). At larger 
values of k, the DS is always subject to an oscillatory 
instability similar to that found for regular dark optical 
solitons [23]. Nevertheless, this instability is character- 
ized by a small growth rate, whose maximum value is 
UJ^ «0.04 (at fc« 0.024). 

An example of the weakly unstable DS for k =- 0.001 



can be seen in the middle panels of Fig. 1, and an exam- 
ple close to the termination of the branch (at k — 0.119) 
is shown in the bottom panels. Note that the instabil- 
ity is manifested by the presence of a nonzero imaginary 
part of the eigenfrequencies. 




FIG. 2. The top panels in Fig. 2 are equivalent to those 
in Fig. 1, but now the strength of the trapping potential is 
fixed, k — 0.001, while the coupling constant C is varied. The 
middle and bottom panels show, respectively, examples of the 
spatial profiles and eigenfrequency spectral planes of unstable 
(for C = 0.25) and stable (for C = 0.01) DSs. 

Next, we consider the case where the strength of the 
parabolic trap is kept fixed, k — 0.001, while the coupling 
constant C is varied in the interval < C < 1. In this 
case, the DS configuration on top of the TF background 
has been obtained for every value of C, see Fig. 2, down 
to C = corresponding to the so-called anti-continuum 
(AC) limit. The norm of the solution (top left panel of 
Fig. 2) is a slowly (almost linearly) decreasing function 
of C, which can be easily understood. Indeed, since the 
effective size of the TF distribution is kept constant (the 
trap strength k is fixed), the DS placed at the center of 
the condensate expands to a larger number of lattice sites 
as C increases, hence a bigger internal part of the BEC 
cloud is efi'ectively devoid of atoms. The dependence of 
the largest unstable eigenfrequency tOi on C, shown in 
the top right panel of Fig. 2, reveals that there exist a 
critical value, Ccr ~ 0.076, at which a Hamiltonian-Hopf 
bifurcation occurs, as two pairs of eigenvalues with op- 
posite Krein signatures [24] collide and bifurcate into a 
complex quartet. 

Notice the proximity of the critical value observed 
herein and the one found in Ref. [23] for regular opti- 
cal DSs (without the parabolic trap). This observation 
suggests that the presence of the parabolic trap does not 
significantly affect the threshold of the oscillatory insta- 
bility. 

It is noteworthy that, at larger values of C (see the 
top right panel of Fig. 2), there are windows on the C 
axis (such as e.g., 0.895 < C < 0.935 for the 200-site 
lattice) where the DS is stable. Examples of unstable 



(for C = 0.25) and stable DSs (for C = 0.01), as well 
as the respective spectral planes (w^, w^) of the eigenfre- 
quencies, are shown in the middle and bottom panels of 
Fig. 2. As noted in Ref. [23], the stability windows result 
from the finitencss of the lattice. 

As the DNLS of Eq. (3) is only a tight-binding approx- 
imation of the full continuum model, i.e., the GP equa- 
tion (1), it is necessary to investigate whether and how 
the weak instability of the DSs, demonstrated above in 
the discrete limiting case, manifests itself in the full GP 
equation. To this end, we first consider the case where 
the OL potential is fixed (an analog of the case where C 
is fixed in the discrete model), with strength Vq = 1.5 
and wavelength A = 5. 




FIG. 3. The top panels show the number of atoms rescaled 
(for convenience) by a factor of 10^^''' (left) and the maxi- 
mum of the instability growth rate LJi (right) for the dark soli- 
ton versus the parabolic trap strength fc, as found from the 
continuum Gross-Pitaevskii equation (1) with Vb = 1.5 and 
A = 5. The middle and bottom panels display the spatial pro- 
files for the dark solitons with k = 0.0024 and k = 0.01 (left) 
and the respective stability spectral planes (right), [k = 0.01 
being very close to the termination point of the branch.] 

Similar to what we find for the discrete case, the DS 
on the TF background exists, in the continuum GP equa- 
tion, only up to a critical value of the trap strength, 
fccr ~ 0.01. This is demonstrated in the top left panel 
of Fig. 3, where the number of atoms, which is now de- 
fined as iV = ||w||2 = /_ \'^{^)\ dx , is shown versus k. 
Furthermore, the eigenvalue computation shows that this 
family of continuum DSs is always unstable (although in 
other cases DSs in the continuum model may be stable, 
see below) . However the instability is extremely weak, as 
its largest growth rate max(ct;i) « 0.015 (at k k, 0.024). 
The shape of typical continuum DSs are shown in the 
middle and bottom panels of Fig. 3 together with the 
corresponding spectral planes. 

Next, we consider the case when the strength of the 
magnetic trap, in the continuum GP equation, is fixed 
at fc = 0.001, while the wavelength A of the OL is varied 



(the strength of the OL is again Vb — 1.5). In this case, 
the DS on the TF background exists for every value of 
A in an interval 3 < A < 5, which is chosen to display 
the present case. In particular, the top left panel in Fig. 
4 shows the number of atoms in the DS as a function of 
A. Furthermore, the dependence of the largest instabil- 
ity growth rate uji on A (top right panel in Fig. 4) reveals 
that a Hamiltonian-Hopf bifurcation occurs at a critical 
point, Acr ~ 4, so that the DSs are stable for A < Acr 
and unstable otherwise. Examples of stable and unstable 
DSs for A = 3 and A = 5 are shown in the middle and 
bottom panels in Fig. 4, together with their respective 
spectral planes. We have verified that for A < 3 the DS 
retains its stability, while for A > 5 it remains unstable 
up to A < Acr ~ 5.45 and restabilizes for larger values of 
A. After running simulations at many other parameter 
values, we have concluded that the OL periodicity A is 
a crucial factor in determining the stability of the DS in 
the continuum GP equation, while the strength of the 
parabolic trap k plays a lesser role. 
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FIG. 4. The top panels show the number of atoms in the 
continuum dark soliton (rescaled by the factor IQ^^''^; left) 
and the largest instability growth rate uii (right) versus the 
OL wavelength A for the GP equation with Vo = 1.5 and 
k — 0.001. The middle and bottom panels show the spatial 
profile of the dark solitons for A = 3 and A = 5 (left) and 
their respective spectral planes (right). 

Finally, an additional factor of relevance in determin- 
ing stability is the amplitude Vb of the optical lattice. We 
have found (data not shown) that variation of Vb results 
in the presence of windows of stability for a fixed value 
of A. For instance for A = 5, the soliton is stable in the 
absence of the optical lattice (i.e., for Vq = 0) as well as 
e.g., in the intervals Vb G [0.04,0.13] or Vb £ [1.17, 1.38], 
while it is unstable for intermediate values. 

In the case when solitons are unstable, it is desirable 
to directly simulate the full nonlinear equation, in order 
to observe the evolution of the soliton under the effect 
of the instability. We find that the instability sets in at 
large times, on the order of O(IOO), and manifests itself 



as follows: the DS, which is initially placed at the bottom 
of the composite potential (magnetic trap and OL), per- 
forms small-amplitude oscillations inside the respective 
small well of the OL potential. Then, as time passes, the 
soliton starts to emit small amounts of radiation, and, as 
a result, it gains enough kinetic energy to perform larger- 
amplitude quasi-periodic oscillations around the center of 
the condensate. This behavior is demonstrated in the left 
panel of Fig. 5, where two snapshots of the density profile 
of the DS are shown at i = and t — 700 for the unstable 
DS of the bottom panel of Fig. 4. To trigger the onset of 
instability, a uniformly distributed random perturbation 
of amplitude 0.01 is added to the initial configuration. It 
can be clearly observed that, in the final state, the DS 
has shifted its position from the center of the condensate 
(x = 0). Quasi-periodic oscillations of the density at the 
point a; = 0, induced by the instability, arc shown in the 
right panel of Fig. 5. 
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FIG. 5. Left: The profile of the dark soliton at i = 
(dashed line) and t = 700 (solid line). Right: Time evolu- 
tion of the density, |w(a;,t)| , 
Vo = 1.5, A = 5 and fc = 0.001. 



at the central point, a; = 0. 



In conclusion, we investigate the existence and stabil- 
ity of stationary dark solitons in repulsive BECs, trapped 
in an optical lattice (OL). The consideration is performed 
in the framework of both the continuous Gross-Pitaevskii 
(GP) equation and its discrete tight-binding dynamical- 
lattice counterpart, taking into regard the effect of the 
external magnetic trap and OL. We find that in the dis- 
crete model the dark solitons are, generally, subject to 
a weak oscillatory instability. In the full GP equation, 
dark solitons may be stable and the stability is chiefly 
determined by the period and amplitude of the OL. If 
the oscillatory instability is present, it sets in at large 
times, which attests to the robustness of dark solitons. 
The instability eventually manifests itself as a shift of the 
dark soliton from the center of the condensate, which is 
accompanied by quasi-periodic oscillations. This work 
can be naturally extended for two-dimensional vortices. 
First steps in that direction were made in Refs. [25] and 
[26]. 
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